
# ------------------------------------------------------------------------------
# Replication Materials
# 
# title: Coalition Size, Direct Democracy, and Public Spending
# journal: Journal of Public Policy
# authors: Patrick Emmenegger, Lucas Leemann, and André Walter
# date: Sept 2020
# ------------------------------------------------------------------------------




#### Result Stability




pdf("out/Time_stable_test.pdf",width=19,height=18)
par(family="CMU Serif", mar=c(5.1, 5.1, 0.5, 2.1), mfrow=c(3,1), cex=2)


move.d <- 0.4
plot(0,0,xlim=c(1870,2005), ylim=c(-1,1), bty="n", xlab="",ylab=expression(hat(beta)[early]-hat(beta)[late]))
N.sim <- 2000
for (y in seq(1870,2000,by = 10)){
  
  findat_time1 <- findat_oLG[findat_oLG$year<=y,]
  findat_time2 <- findat_oLG[findat_oLG$year>y,]
  
  m1 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time1)
  m2 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time2)
  
  
  b1 <- coef(m1)[c("numpar:index_law_init","index_law_init")]
  b2 <- coef(m2)[c("numpar:index_law_init","index_law_init")]
  V1 <- vcov(m1)[c("numpar:index_law_init","index_law_init"),c("numpar:index_law_init","index_law_init")]
  V2 <- vcov(m2)[c("numpar:index_law_init","index_law_init"),c("numpar:index_law_init","index_law_init")]
  
  
  BETA1 <- mvrnorm(N.sim, b1, V1)
  BETA2 <- mvrnorm(N.sim, b2, V2)
  
  B1 <- BETA1[,c("numpar:index_law_init","index_law_init")]
  BB1 <- B1[,2] + 5* B1[,1]
  
  B2 <- BETA2[,c("numpar:index_law_init","index_law_init")]
  BB2 <- B2[,2] + 5* B2[,1]
  
  
  dd <- (BB1) - (BB2)
  
  # points(rep(y,N.sim)+rnorm(N.sim,sd=.25),
  #        dd, 
  #        col=rgb(100,0,0,4,maxColorValue = 255), pch=19)
  abline(h=0, col=rgb(100,0,0,255,maxColorValue = 255), lwd=2)
  
  QQ <-  quantile(dd,c(0.005,0.025,0.975,0.995))
  print(QQ)
  segments(y-move.d,QQ[2],y-move.d,QQ[3],lwd=3,col=rgb(100,0,0,255,maxColorValue = 255))
  segments(y-move.d,QQ[1],y-move.d,QQ[4],lwd=1,col=rgb(100,0,0,255,maxColorValue = 255))
  
}






move.d <- 0.4
plot(0,0,xlim=c(1870,2005), ylim=c(-.5,.5), bty="n", xlab="",ylab=expression(hat(beta)[early]-hat(beta)[late]))
N.sim <- 2000
for (y in seq(1870,2000,by = 10)){
  
  findat_time1 <- findat_oLG[findat_oLG$year<=y,]
  findat_time2 <- findat_oLG[findat_oLG$year>y,]
  
  m1 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time1)
  m2 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time2)
  
  
  b1 <- coef(m1)[c("numpar:index_law_ref","index_law_ref")]
  b2 <- coef(m2)[c("numpar:index_law_ref","index_law_ref")]
  V1 <- vcov(m1)[c("numpar:index_law_ref","index_law_ref"),c("numpar:index_law_ref","index_law_ref")]
  V2 <- vcov(m2)[c("numpar:index_law_ref","index_law_ref"),c("numpar:index_law_ref","index_law_ref")]
  
  
  BETA1 <- mvrnorm(N.sim, b1, V1)
  BETA2 <- mvrnorm(N.sim, b2, V2)
  
  B1 <- BETA1[,c("numpar:index_law_ref","index_law_ref")]
  BB1 <- B1[,2] + 5* B1[,1]
  
  B2 <- BETA2[,c("numpar:index_law_ref","index_law_ref")]
  BB2 <- B2[,2] + 5* B2[,1]
  
  dd <- (BB1) - (BB2)
  
  # points(rep(y,N.sim)+rnorm(N.sim,sd=.25),
  #        dd, 
  #        col=rgb(100,0,0,4,maxColorValue = 255), pch=19)
  abline(h=0, col=rgb(100,0,0,255,maxColorValue = 255), lwd=2)
  
  QQ <-  quantile(dd,c(0.005,0.025,0.975,0.995))
  print(QQ)
  segments(y-move.d,QQ[2],y-move.d,QQ[3],lwd=3,col=rgb(100,0,0,255,maxColorValue = 255))
  segments(y-move.d,QQ[1],y-move.d,QQ[4],lwd=1,col=rgb(100,0,0,255,maxColorValue = 255))
  
}


move.d <- 0.4
plot(0,0,xlim=c(1870,2005), ylim=c(-10,10), bty="n", xlab="",ylab=expression(hat(beta)[early]-hat(beta)[late]))
N.sim <- 2000
for (y in seq(1870,2000,by = 10)){
  
  findat_time1 <- findat_oLG[findat_oLG$year<=y,]
  findat_time2 <- findat_oLG[findat_oLG$year>y,]
  
  m1 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time1)
  m2 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time2)
  
  
  b1 <- coef(m1)[c("numpar:refmon_adj_imp0","refmon_adj_imp0")]
  b2 <- coef(m2)[c("numpar:refmon_adj_imp0","refmon_adj_imp0")]
  V1 <- vcov(m1)[c("numpar:refmon_adj_imp0","refmon_adj_imp0"),c("numpar:refmon_adj_imp0","refmon_adj_imp0")]
  V2 <- vcov(m2)[c("numpar:refmon_adj_imp0","refmon_adj_imp0"),c("numpar:refmon_adj_imp0","refmon_adj_imp0")]
  
  
  BETA1 <- mvrnorm(N.sim, b1, V1)
  BETA2 <- mvrnorm(N.sim, b2, V2)
  
  B1 <- BETA1[,c("numpar:refmon_adj_imp0","refmon_adj_imp0")]
  BB1 <- B1[,2] + 5* B1[,1]
  
  B2 <- BETA2[,c("numpar:refmon_adj_imp0","refmon_adj_imp0")]
  BB2 <- B2[,2] + 5* B2[,1]
  
  
  dd <- (BB1) - (BB2)
  
  # points(rep(y,N.sim)+rnorm(N.sim,sd=.25),
  #        dd, 
  #        col=rgb(100,0,0,4,maxColorValue = 255), pch=19)
  abline(h=0, col=rgb(100,0,0,255,maxColorValue = 255), lwd=2)
  
  QQ <-  quantile(dd,c(0.005,0.025,0.975,0.995))
  print(QQ)
  segments(y-move.d,QQ[2],y-move.d,QQ[3],lwd=3,col=rgb(100,0,0,255,maxColorValue = 255))
  segments(y-move.d,QQ[1],y-move.d,QQ[4],lwd=1,col=rgb(100,0,0,255,maxColorValue = 255))
  
}



dev.off()








pdf("out/Time_stable_woWWII.pdf",width=19,height=6)
par(family="CMU Serif", mar=c(5.1, 5.1, 0, 2.1))
move.d <- 0.4
plot(0,0,xlim=c(1870,2005), ylim=c(-1,1), bty="n", xlab="",ylab=expression(hat(beta)[early]-hat(beta)[late]))
N.sim <- 2000
for (y in seq(1870,2000,by = 10)){
  
  findat_time1 <- findat_oLG[findat_oLG$year<=y,]
  findat_time2 <- findat_oLG[findat_oLG$year>y,]
  
  m1 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time1[findat_time1$year!=1935 &findat_time1$year!=1940 & findat_time1$year!=1945 ,])
  m2 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix, findat_time2[findat_time2$year!=1935 &findat_time2$year!=1940 & findat_time2$year!=1945 ,])
  
  
  b1 <- coef(m1)[c("numpar:index_law_init","index_law_init")]
  b2 <- coef(m2)[c("numpar:index_law_init","index_law_init")]
  V1 <- vcov(m1)[c("numpar:index_law_init","index_law_init"),c("numpar:index_law_init","index_law_init")]
  V2 <- vcov(m2)[c("numpar:index_law_init","index_law_init"),c("numpar:index_law_init","index_law_init")]
  
  
  BETA1 <- mvrnorm(N.sim, b1, V1)
  BETA2 <- mvrnorm(N.sim, b2, V2)
  
  B1 <- BETA1[,c("numpar:index_law_init","index_law_init")]
  BB1 <- B1[,2] + 5* B1[,1]
  
  B2 <- BETA2[,c("numpar:index_law_init","index_law_init")]
  BB2 <- B2[,2] + 5* B2[,1]
  
  
  dd <- (BB1) - (BB2)
  
  # points(rep(y,N.sim)+rnorm(N.sim,sd=.25),
  #        dd, 
  #        col=rgb(100,0,0,4,maxColorValue = 255), pch=19)
  abline(h=0, col=rgb(100,0,0,255,maxColorValue = 255), lwd=2)
  
  QQ <-  quantile(dd,c(0.005,0.025,0.975,0.995))
  print(QQ)
  segments(y-move.d,QQ[2],y-move.d,QQ[3],lwd=3,col=rgb(100,0,0,255,maxColorValue = 255))
  segments(y-move.d,QQ[1],y-move.d,QQ[4],lwd=1,col=rgb(100,0,0,255,maxColorValue = 255))
  print(y)
}

dev.off()









outi <- round(psych::describe(findat_oLG[,c("dlnexp_pc","year","numpar", "llnexp_pc", "index_law_init","index_law_ref",
                                            "refmon_adj_imp0", "sec_sh", "fir_sh", "deprat", "babydeathratio_impute", 
                                            "left", "PR", "phydens", "lnpop", "lnfed_mix", "fisrul")], skew=FALSE, check=TRUE),1)

as.matrix(outi)
print(xtable(outi), file="out/table_3.tex")



mm1 <- brm(dlnexp_pc~factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix+(1|canton) + (-1+index_law_init:numpar+index_law_init|year) , 
           data = findat_oLG, chains = 4, cores= 2)

mm2 <- brm(dlnexp_pc~factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix+(1|canton) + (-1+index_law_ref:numpar+index_law_ref|year) , 
           data = findat_oLG, chains = 4, cores= 2)

mm3 <- brm(dlnexp_pc~factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
             fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix+(1|canton) + (-1+refmon_adj_imp0:numpar+refmon_adj_imp0|year) , 
           data = findat_oLG, chains = 4, cores= 2)

pdf("out/Time_stable_RE_test.pdf",width=22,height=22)
par(family="CMU Serif", mar=c(5.1, 5.1, 0.5, 2.1), mfrow=c(3,1), cex=2)
plot(seq(1860,2015,by = 5), ranef(mm1)$year[,1,1], bty="n",
     ylim=c(-0.06,0.06), ylab="Random Effect Part of the Coefficient",
     xlab="", xaxt="n", pch=19, col=rgb(100,0,0,255,maxColorValue = 255))
segments(seq(1860,2015,by = 5),
         ranef(mm1)$year[,3,1],
         seq(1860,2015,by = 5),
         ranef(mm1)$year[,4,1], 
         col=rgb(100,0,0,255,maxColorValue = 255), lwd=1.5)
abline(h=0, lwd=3)
axis(1,seq(1860,2010,by = 10),labels=seq(1860,2010,by = 10))



#pdf("Time_stable_RE_LR.pdf",width=19,height=8)
#par(family="CMU Serif", mar=c(5.1, 5.1, 0, 2.1))
plot(seq(1860,2015,by = 5), ranef(mm2)$year[,1,1], bty="n",
     ylim=c(-0.06,0.06), ylab="Random Effect Part of the Coefficient",
     xlab="", xaxt="n", pch=19, col=rgb(100,0,0,255,maxColorValue = 255))
segments(seq(1860,2015,by = 5),
         ranef(mm2)$year[,3,1],
         seq(1860,2015,by = 5),
         ranef(mm2)$year[,4,1], 
         col=rgb(100,0,0,255,maxColorValue = 255), lwd=1.5)
abline(h=0, lwd=3)
axis(1,seq(1860,2010,by = 10),labels=seq(1860,2010,by = 10))
#dev.off()




#pdf("Time_stable_RE_FR.pdf",width=19,height=8)
#par(family="CMU Serif", mar=c(5.1, 5.1, 0, 2.1))
plot(seq(1860,2015,by = 5), ranef(mm3)$year[,1,1], bty="n",
     ylim=c(-2,2), ylab="Random Effect Part of the Coefficient",
     xlab="", xaxt="n", pch=19, col=rgb(100,0,0,255,maxColorValue = 255))
segments(seq(1860,2015,by = 5),
         ranef(mm3)$year[,3,1],
         seq(1860,2015,by = 5),
         ranef(mm3)$year[,4,1], 
         col=rgb(100,0,0,255,maxColorValue = 255), lwd=1.5)
abline(h=0, lwd=3)
axis(1,seq(1860,2010,by = 10),labels=seq(1860,2010,by = 10))
dev.off()



